﻿
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;

namespace LxCFD
{
    public class PressureCorrectionBasic : EquationBase
    {
        public Vector3d[] Dc;//Dc，见591。用以求解Df,DDf。DDf见P596//Dc在保存的时候，剔除松弛的影响！

        //not local
        public VelocityVector VelocityVector { get; set; }
        protected double[] MassFlux;
        protected double[] MassImbalance;


        public PressureSIMPLE PressureSIMPLE;
        protected double[] Pressure;

        protected Value1Base<double> Density;
        public SIMPLEMethods SIMPLEMethod { get; set; }

        #region 公共
        public override void CheckCreate(FvZone zone)
        {
            base.CheckCreate(zone);
            LxTool.CheckCreateArray(ref Dc, Cells.Length);
        }
        public void SetZone(double pressureRelaxFactor, SIMPLEMethods solvingMethod, VelocityVector velocityVector, PressureSIMPLE pressureInSimple, Value1Base<double> density)
        {
            clear_eFlags();
            this.RelaxFactor = pressureRelaxFactor;
            this.SIMPLEMethod = solvingMethod;

            this.VelocityVector = velocityVector;
            this.MassFlux = velocityVector.MassFlux;
            this.MassImbalance = velocityVector.MassImbalance;
 

            this.PressureSIMPLE = pressureInSimple;
            this.Pressure = pressureInSimple.x;

            this.Density = density;
        }

        public void SetBound_AfterSettingVelocityVector()
        {
            foreach (var b in VelocityVector.NoSlipWalls)
            {
                e_ALUO0.Add(b.ThreadID);
            }
            foreach (var b in VelocityVector.Symmetries)
            {
                e_ALUO0.Add(b.ThreadID);
            }
            foreach (var b in VelocityVector.VelocityInlets)
            {
                e_ALUO0.Add(b.ThreadID);
            }
            foreach (var b in VelocityVector.PressureInOuts)
            {
                e_FixValue.Add(b.ThreadID);
            }
        }
        #endregion
        #region 外部调用
        public override void Set_x_Bound(int taskIndex)//压力修正初始值都为0
        {
            TotalIndex.CellIDs.GetSeg(out int i1, out int endp1, taskIndex);
            for (; i1 < endp1; i1++)
            {
                x[i1] = 0;
            }
        }
        public override void SetAb_Step1(int taskIndex)
        {
            set_ALU_inner(taskIndex);
            foreach (var b in VelocityVector.NoSlipWalls)
            {
                set_ALU_0(b.ThreadIndex, taskIndex);
            }
            foreach (var b in VelocityVector.Symmetries)
            {
                set_ALU_0(b.ThreadIndex, taskIndex);
            }
            foreach (var b in VelocityVector.VelocityInlets)
            {
                set_ALU_0(b.ThreadIndex, taskIndex);
            }
            foreach (var b in VelocityVector.PressureInOuts)
            {
                set_ALU_bound(b.ThreadIndex, taskIndex);
            }
        }
        public override void SetAb_Step2(int taskIndex)
        {
            InnerIndex.CellIDs.GetSeg(out int i1, out int endp1, taskIndex);
            for (; i1 < endp1; i1++)
            {
                b[i1] = -MassImbalance[i1];
                ADiag[i1] = -Cells[i1].CalAllAdd(ALUN, ALUO);
            }
        }
        public void SetAb_PCor0Point(int CellIndex)
        {
            ref FvCell cell = ref Cells[CellIndex];
            foreach (var f in cell.OFaces)
            {
                ALUO[f] = 0;
            }
            foreach (var f in cell.NFaces)
            {
                ALUN[f] = 0;
            }
            b[CellIndex] = 0;
        }
        public override void Relax(int taskIndex)
        {
        }
        public override void Update_x_Bound(int taskIndex)
        {
            foreach (var b in VelocityVector.NoSlipWalls)
            {
                predict_Bound(b.ThreadIndex, taskIndex);
            }
            foreach (var b in VelocityVector.Symmetries)
            {
                predict_Bound(b.ThreadIndex, taskIndex);
            }
            foreach (var b in VelocityVector.VelocityInlets)
            {
                predict_Bound(b.ThreadIndex, taskIndex);
            }
        }
        public void ApplyPCorToUVWP(int taskIndex)
        {
            addCorrectionToUVW(InnerIndex, taskIndex);
            addCorrectionToPressure(InnerIndex, taskIndex);
            addPressureCorrectionToMassFlux_inner(InnerIndex, taskIndex);
            foreach (var b in VelocityVector.PressureInOuts)
            {
                addPressureCorrectionToMassFlux_bound(b.ThreadIndex, taskIndex);
            }
        }
        public void ApplyPCorToUVWP_NoRelax(int taskIndex)
        {
            addCorrectionToUVW(InnerIndex, taskIndex);
            addCorrectionToPressure_NoRelax(InnerIndex, taskIndex);
            addPressureCorrectionToMassFlux_inner(InnerIndex, taskIndex);
            foreach (var b in VelocityVector.PressureInOuts)
            {
                addPressureCorrectionToMassFlux_bound(b.ThreadIndex, taskIndex);
            }
        }
        public void CalDc_ADiagNotRelaxed(int taskIndex)//如果U V W的系数矩阵未进行过松弛//添加源项可能改动A_Diag，所以计算完A_Diag后立即计算Dc
        {
            var U_A_Diag = VelocityVector.U.ADiag;
            var V_A_Diag = VelocityVector.V.ADiag;
            var W_A_Diag = VelocityVector.W.ADiag;
            InnerIndex.CellIDs.GetSeg(out int i1, out int endp1, taskIndex);
            for (; i1 < endp1; i1++)
            {
                double volumn = Cells[i1].Volumn;
                Dc[i1] = new Vector3d(volumn / U_A_Diag[i1], volumn / V_A_Diag[i1], volumn / W_A_Diag[i1]);
            }
        }
        public virtual void Setb_Step1_Cor(int taskIndex)
        {
            throw new NotImplementedException();
        }
        public virtual void Setb_Step2_Cor(int taskIndex)
        {
            throw new NotImplementedException();
        }
        #endregion
        protected virtual void set_ALU_inner(int taskIndex)
        {
            InnerIndex.FaceIDs.GetSeg(out int i1, out int endp1, taskIndex);
            var mface = new PhysicalFaceOverRelax();
            for (; i1 < endp1; i1++)
            {
                ref FvFace face = ref Faces[i1];
                Vector3d Df_average = face.CalAverage(Dc);
                mface.SetFace(ref face, ref Df_average);
                ALUO[i1] = ALUN[i1] = -Density.GetFaceValue(i1) * mface.gDiff;
            }
        }
        protected virtual void set_ALU_bound(IndexNumber indexNumber, int taskIndex)
        {
            indexNumber.FaceIDs.GetSeg(out int i1, out int endp1, taskIndex);
            var mface = new PhysicalFaceOverRelax();
            for (; i1 < endp1; i1++)
            {
                ref FvFace face = ref Faces[i1];
                Vector3d Df_average = Dc[face.OSideCell];
                mface.SetFace(ref face, ref Df_average);
                ALUO[i1] = ALUN[i1] = -Density.GetFaceValue(i1) * mface.gDiff;
            }
        }
        protected void predict_Bound(IndexNumber indexNumber, int taskIndex)
        {
            indexNumber.FaceIDs.GetSeg(out int i1, out int endp1, taskIndex);
            for (; i1 < endp1; i1++)
            {
                x[Faces[i1].NSideCell] = x[Faces[i1].OSideCell];
            }
        }
        #region 进行修正
        protected void addCorrectionToUVW(IndexNumber indexNumber, int taskIndex)
        {
            var U_x = VelocityVector.U.x;
            var V_x = VelocityVector.V.x;
            var W_x = VelocityVector.W.x;
            indexNumber.CellIDs.GetSeg(out int i1, out int endp1, taskIndex);
            for (; i1 < endp1; i1++)
            {
                Vector3d velocityCorrection = Vector3d.XXYYZZProduct(Dc[i1], Grad[i1]);
                U_x[i1] -= velocityCorrection.X;
                V_x[i1] -= velocityCorrection.Y;
                W_x[i1] -= velocityCorrection.Z;
            }
        }
        protected void addCorrectionToPressure(IndexNumber indexNumber, int taskIndex)
        {
            indexNumber.CellIDs.GetSeg(out int i1, out int endp1, taskIndex);
            for (; i1 < endp1; i1++)
            {
                Pressure[i1] += RelaxFactor * x[i1];
            }
        }
        protected void addCorrectionToPressure_NoRelax(IndexNumber indexNumber, int taskIndex)
        {
            indexNumber.CellIDs.GetSeg(out int i1, out int endp1, taskIndex);
            for (; i1 < endp1; i1++)
            {
                Pressure[i1] += x[i1];
            }
        }
        protected virtual void addPressureCorrectionToMassFlux_inner(IndexNumber indexNumber, int taskIndex)//需要对公式校核以下，这个公式理论上应与RhieChow一致
        {
            indexNumber.FaceIDs.GetSeg(out int i1, out int endp1, taskIndex);
            for (; i1 < endp1; i1++)
            {
                ref FvFace face = ref Faces[i1];
                Vector3d Df_average = face.CalAverage(Dc);
                Vector3d gradPf = face.CalCFCorrectedGradient_inner(Grad, x);
                double vf = Vector3d.DotProduct(Vector3d.XXYYZZProduct(Df_average, gradPf), face.Normal);//P596
                MassFlux[i1] -= Density.GetFaceValue(i1) * vf;
            }
        }
        protected virtual void addPressureCorrectionToMassFlux_bound(IndexNumber indexNumber, int taskIndex)//需要对公式校核以下，这个公式理论上应与RhieChow一致
        {
            indexNumber.FaceIDs.GetSeg(out int i1, out int endp1, taskIndex);
            for (; i1 < endp1; i1++)
            {
                ref FvFace face = ref Faces[i1];
                Vector3d Df_average = Dc[face.OSideCell];
                Vector3d gradPf = face.CalCFCorrectedGradient_bound(Grad, x);
                double vf = Vector3d.DotProduct(Vector3d.XXYYZZProduct(Df_average, gradPf), face.Normal);//P620,公式是否正确？？？
                MassFlux[i1] -= Density.GetFaceValue(i1) * vf;
            }
        }
        #endregion
    }
}
//protected double[] U_x;
//protected double[] V_x;
//protected double[] W_x;
//protected double[] U_A_Diag;
//protected double[] V_A_Diag;
//protected double[] W_A_Diag;
//this.U_x = velocityVector.U.x;
//this.V_x = velocityVector.V.x;
//this.W_x = velocityVector.W.x;
//this.U_A_Diag = velocityVector.U.ADiag;
//this.V_A_Diag = velocityVector.V.ADiag;
//this.W_A_Diag = velocityVector.W.ADiag;
//public List<BoundaryHeadBase> FixedPressure { get; set; } = new List<BoundaryHeadBase>();
//public List<BoundaryHeadBase> ChangablePressure { get; set; } = new List<BoundaryHeadBase>();
//FixedPressure.Clear();
//ChangablePressure.Clear();
//public void SetBound_FixP(FvThread thread)
//{
//    BoundaryHeadBase b = new BoundaryHeadBase();
//    b.SetHead(thread);

//    FixedPressure.Add(b);
//    e_FixValue.Add(thread.ID);
//}
//public void SetBound_ChangableP(FvThread thread)
//{
//    BoundaryHeadBase b = new BoundaryHeadBase();
//    b.SetHead(thread);

//    ChangablePressure.Add(b);
//    e_ALUO0.Add(thread.ID);
//}
//addCorrectionToUVWP(InnerIndex, taskIndex);
//addCorrectionToUVWP_NoRelax(InnerIndex, taskIndex);
//protected void addCorrectionToUVWP_NoRelax(IndexNumber indexNumber, int taskIndex)
//{
//    indexNumber.CellIDs.GetSeg(out int i1, out int endp1, taskIndex);
//    for (; i1 < endp1; i1++)
//    {
//        //UVW
//        Vector3d velocityCorrection = Vector3d.XXYYZZProduct(Dc[i1], Grad[i1]);
//        U_x[i1] -= velocityCorrection.X;
//        V_x[i1] -= velocityCorrection.Y;
//        W_x[i1] -= velocityCorrection.Z;
//        //P
//        Pressure[i1] += x[i1];
//    }
//}
//protected void addCorrectionToUVWP(IndexNumber indexNumber, int taskIndex)
//{
//    indexNumber.CellIDs.GetSeg(out int i1, out int endp1, taskIndex);
//    for (; i1 < endp1; i1++)
//    {
//        //UVW
//        Vector3d velocityCorrection = Vector3d.XXYYZZProduct(Dc[i1], Grad[i1]);
//        U_x[i1] -= velocityCorrection.X;
//        V_x[i1] -= velocityCorrection.Y;
//        W_x[i1] -= velocityCorrection.Z;
//        //P
//        Pressure[i1] += RelaxFactor * x[i1];
//    }
//}